# Loads data
load("bc_bear_occurrences.Rda")
load("BC_Covariates.Rda")
# Generates summary of loaded data
summary(DATA)
## Length Class Mode
## Window 1 SpatialPolygons S4
## Elevation 10 im list
## Forest 10 im list
## HFI 10 im list
## Dist_Water 10 im list
# Extracts location columns
bears_loc <- occ_data[, c("decimalLongitude", "decimalLatitude", "month", "year")]
bears_loc_filtered <- subset(bears_loc, year %in% c(2020, 2021, 2022, 2023, 2024))
# Creates sf object
bears_sf <- st_as_sf(
bears_loc_filtered,
coords = c("decimalLongitude", "decimalLatitude"),
crs = 4326 # WGS84 (longitude/latitude)
)
# BC Albers projection
bears_sf_proj <- st_transform(bears_sf, crs = 3005)
# Extracts coordinates
coords <- st_coordinates(bears_sf_proj)
# Extracts BC window
window_sf <- st_as_sf(DATA$Window)
window_proj <- st_transform(window_sf, crs = 3005)
window <- as.owin(window_proj)
# Creates ppp object
bears_ppp <- ppp(
x = coords[,1],
y = coords[,2],
window = window
)
## Warning: 187 points were rejected as lying outside the specified window
## Warning: data contain duplicated points
# Plots bear occurrences
plot(bears_ppp,
pch = 21,
main = "Black bear occurrences in BC, 2024")
## Warning in plot.ppp(bears_ppp, pch = 21, main = "Black bear occurrences in BC,
## 2024"): 187 illegal points also plotted
# Defines the seasons and corresponding colors
seasons <- list(
winter = c(12, 1, 2),
spring = c(3, 4, 5),
summer = c(6, 7, 8),
autumn = c(9, 10, 11)
)
# Creates an empty list to store the ppp objects for each season
ppp_list <- list()
# Creates an empty list to store the filtered data for each season
bears_sf_list <- list()
# Loops over seasons
for (i in 1:length(seasons)) {
# Filters for each season
season_name <- names(seasons)[i]
season_months <- seasons[[i]]
# Filters bears_loc_filtered by the season's months
bears_loc_season <- bears_loc_filtered[bears_loc_filtered$month %in% season_months, ]
# Prints the number of observations for this season
cat(season_name, ": ", nrow(bears_loc_season), " observations\n", sep="")
# Creates sf object for the filtered season
bears_sf_season <- st_as_sf(
bears_loc_season,
coords = c("decimalLongitude", "decimalLatitude"),
crs = 4326 # WGS84 (longitude/latitude)
)
# Transforms to BC Albers projection
bears_sf_season_proj <- st_transform(bears_sf_season, crs = 3005)
# Stores sf object in the list
bears_sf_list[[season_name]] <- bears_sf_season_proj
# Extracts x, y coordinates
coords <- st_coordinates(bears_sf_season_proj)
# Creates ppp object for each season
suppressWarnings(
bears_ppp_season <- ppp(
x = coords[, 1],
y = coords[, 2],
window = window
)
)
# Stores ppp object in the list
ppp_list[[season_name]] <- bears_ppp_season
}
## winter: 76 observations
## spring: 808 observations
## summer: 1862 observations
## autumn: 841 observations
# Plots four maps in a 2x2 grid
par(mfrow=c(2,2), mar=c(1,1,1,1))
# Loops over ppp objects and plots
for (i in 1:length(ppp_list)) {
season_name <- names(ppp_list)[i]
suppressWarnings(
plot(ppp_list[[season_name]],
main = paste(season_name, "(2024)"),
pch = 21)
)
}
# Visualizes covariates
plot(DATA$Elevation)
plot(DATA$Forest)
plot(DATA$HFI)
plot(DATA$Dist_Water)
# Isolates covariates
elev <- DATA$Elevation
cover <- DATA$Forest
dist_water <- DATA$Dist_Water
hfi <- DATA$HFI
# Intensity as a function of elevation overall
rho <- rhohat(bears_ppp, elev)
plot(rho,
xlim=c(0, max(elev)),
main="Estimated Rho vs. Elevation")
# Intensity as a function of elevation broken up by seasons
par(mfrow=c(2,2), mar=c(2,2,2,2))
for (i in 1:length(ppp_list)) {
season_name <- names(ppp_list)[i]
rho <- rhohat(ppp_list[[season_name]], elev)
plot(rho,
xlim=c(0, max(elev)),
main = paste(season_name))
}
# Loads elevation raster
elev <- DATA$Elevation
# Converts im to raster
im2raster <- function(im_obj) {
m <- as.matrix(im_obj)
r <- raster(m)
ext <- c(im_obj$xrange[1] - im_obj$xstep/2,
im_obj$xrange[2] + im_obj$xstep/2,
im_obj$yrange[1] - im_obj$ystep/2,
im_obj$yrange[2] + im_obj$ystep/2)
extent(r) <- ext
return(r)
}
elev_raster <- im2raster(elev)
# Extracts elevation values at occurrence locations
bears_sf_proj$elev_value <- extract(elev_raster, bears_sf_proj)
summary(bears_sf_proj$elev_value)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 11.68 461.09 644.88 754.84 937.91 2308.18 1499
KDE to estimate the concentration of bear occurrences based on elevation values. The elevation data for bear locations ranges from a minimum of 11.68 meters to a maximum of 2308.18 meters, with most occurrences concentrated between 461.09 meters and 937.91 meters. The mean elevation is 754.84 meters, slightly higher than the median, indicating a skew towards higher elevations.
# KDE for bear locations
kde_b <- density(bears_sf_proj$elev_value, na.rm = TRUE, bw = "SJ-dpi")
# Summarizes
summary(bears_sf_proj$elev_value)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 11.68 461.09 644.88 754.84 937.91 2308.18 1499
plot(kde_b,
main = "Kernel Density Estimate of Elevation at Bear Locations",
xlab = "Elevation (meters)",
ylab = "Density",
col = "#2C6B2F",
lwd = 2)
# NA with mean elevation
mean_elev <- mean(as.vector(as.matrix(elev)), na.rm = TRUE)
elev_clean <- eval.im(ifelse(is.na(elev), mean_elev, elev))
# Establishes elevation model
elevation_model <- ppm(bears_ppp, ~ elev, covariates = list(elev = elev_clean))
print(elevation_model)
## Nonstationary Poisson process
## Fitted to point pattern dataset 'bears_ppp'
##
## Log intensity: ~elev
##
## Fitted trend coefficients:
## (Intercept) elev
## -17.36871576 -0.00253243
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -17.36871576 2.963479e-02 -17.426798885 -17.310632643 ***
## elev -0.00253243 4.234993e-05 -0.002615434 -0.002449426 ***
## Zval
## (Intercept) -586.09208
## elev -59.79774
# Predicts and plots intensity
pred_intensity_elev <- predict(elevation_model)
plot(pred_intensity_elev,
main = "Predicted Intensity Based on Elevation")
points(bears_ppp,
pch = 20,
cex = 0.4,
col = "#B24A2F")
Fitted model: Intercept: -17.37; Elevation coefficient: -0.0025. The negative coefficient for elevation indicates that, as elevation increases, the expected intensity of bear occurrences decreases. This suggests that bears are less likely to be found at higher elevations.
# Establishes elevation residuals
elevation_residuals <- residuals(elevation_model,
type = "pearson")
# Removes any non-finite residuals
elevation_residuals$v[!is.finite(elevation_residuals$v)] <- NA
# Plots elevation residuals
plot(elevation_residuals,
main = "Residuals - Elevation Model",
na.col = "transparent")
# Elevation model
elevation_model <- ppm(bears_ppp, ~ elev, covariates = list(elev = elev_clean))
# Computes partial residuals for elevation
par_res_elev_all <- parres(elevation_model, "elev")
# Plots partial residuals
plot(par_res_elev_all,
legend = FALSE,
lwd = 2,
main = "Partial Residuals Elevation - All Seasons",
xlab = "Elevation (km)",
ylab = "Partial Residuals")
models_elev_seasonal <- list()
# Loops over season in the ppp_list and models elevation
for (season in names(ppp_list)) {
models_elev_seasonal[[season]] <- ppm(ppp_list[[season]], ~ elev, covariates = list(elev = elev_clean))
}
## Warning: glm.fit: algorithm did not converge
# 2x2 plot layout
par(mfrow = c(2, 2))
# Generates and plots partial residuals per season
for (season in names(models_elev_seasonal)) {
model <- models_elev_seasonal[[season]]
parres_season <- parres(model, "elev")
plot(parres_season,
legend = FALSE,
lwd = 2,
main = paste("Partial Residuals - Elevation", season),
xlab = "Elevation (km)",
ylab = "Partial Residuals")
}
# Seasonal elevation models
models_elev_seasonal <- list()
for (season in names(ppp_list)) {
models_elev_seasonal[[season]] <- ppm(ppp_list[[season]], ~ elev, covariates = list(elev = elev_clean))
cat("Model for", season, ":\n")
print(models_elev_seasonal[[season]])
}
## Warning: glm.fit: algorithm did not converge
## Model for winter :
## Nonstationary Poisson process
## Fitted to point pattern dataset 'ppp_list[[season]]'
##
## Log intensity: ~elev
##
## Fitted trend coefficients:
## (Intercept) elev
## -20.281500393 -0.004896402
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -20.281500393 0.1590970890 -20.593324958 -19.96967583 ***
## elev -0.004896402 0.0004162485 -0.005712234 -0.00408057 ***
## Zval
## (Intercept) -127.47876
## elev -11.76317
## *** Fitting algorithm for 'glm' did not converge ***
## Model for spring :
## Nonstationary Poisson process
## Fitted to point pattern dataset 'ppp_list[[season]]'
##
## Log intensity: ~elev
##
## Fitted trend coefficients:
## (Intercept) elev
## -18.938423756 -0.002390713
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -18.938423756 6.296883e-02 -19.061840399 -18.815007114 ***
## elev -0.002390713 8.708531e-05 -0.002561397 -0.002220029 ***
## Zval
## (Intercept) -300.75870
## elev -27.45254
## Model for summer :
## Nonstationary Poisson process
## Fitted to point pattern dataset 'ppp_list[[season]]'
##
## Log intensity: ~elev
##
## Fitted trend coefficients:
## (Intercept) elev
## -18.213839873 -0.002227723
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -18.213839873 4.265390e-02 -18.297439975 -18.13023977 ***
## elev -0.002227723 5.699771e-05 -0.002339437 -0.00211601 ***
## Zval
## (Intercept) -427.01468
## elev -39.08443
## Model for autumn :
## Nonstationary Poisson process
## Fitted to point pattern dataset 'ppp_list[[season]]'
##
## Log intensity: ~elev
##
## Fitted trend coefficients:
## (Intercept) elev
## -18.514156648 -0.003139775
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -18.514156648 5.671064e-02 -18.625307457 -18.403005838 ***
## elev -0.003139775 9.455484e-05 -0.003325099 -0.002954451 ***
## Zval
## (Intercept) -326.46708
## elev -33.20586
# Combines seasonal data
bears_all_elevation <- do.call(rbind, lapply(names(bears_sf_list), function(season) {
sf_season <- bears_sf_list[[season]]
coords <- st_coordinates(sf_season)
elev_vals <- extract(elev_raster, sf_season)
data.frame(x = coords[, 1], y = coords[, 2],
elev_value = elev_vals, season = season)
}))
bears_all_elevation$season <- as.factor(bears_all_elevation$season)
bears_all_elevation$count <- 1
Winter: Intercept: -20.28; Elevation Coefficient: -0.0049. For winter, the negative coefficient for elevation indicates that bear occurrences decrease as elevation increases. The model did not converge, so caution is needed in interpreting the results.
Spring: Intercept: -18.94; Elevation Coefficient: -0.0024. In spring, the relationship between elevation and bear occurrences is also negative, with a slight decrease in occurrences as elevation increases. However, the model showed some issues with convergence.
Summer: Intercept: -18.21; Elevation Coefficient: -0.0022. Summer shows a similar negative relationship, with bear occurrences decreasing slightly as elevation increases.
Autumn: Intercept: -18.51; Elevation Coefficient: -0.0031. In autumn, bear occurrences also decrease with elevation, with a slightly stronger negative coefficient compared to spring and summer.
All seasons show a negative relationship between bear occurrences and elevation. Winter has the largest negative coefficient (-0.0049), meaning a stronger decrease in occurrences with elevation.
# GAM with elevation and season
gam_elevation <- gam(count ~ s(elev_value) + season, data = bears_all_elevation, family = poisson())
summary(gam_elevation)
##
## Family: poisson
## Link function: log
##
## Formula:
## count ~ s(elev_value) + season
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.729e-17 4.510e-02 0 1
## seasonspring 2.172e-22 6.531e-02 0 1
## seasonsummer 7.375e-22 5.431e-02 0 1
## seasonwinter 8.153e-21 1.544e-01 0 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(elev_value) 1 1.001 0 1
##
## R-sq.(adj) = NaN Deviance explained = NaN%
## UBRE = -0.99521 Scale est. = 1 n = 2088
GAM results suggest that neither elevation nor season has a significant influence on the spatial distribution of bear locations. Warnings about the model fitting process suggest that the results may not be reliable.
# Polynomial ppm model with season
bears_all_elevation <- ppp(
x = bears_all_elevation$x,
y = bears_all_elevation$y,
window = window,
marks = bears_all_elevation$season
)
## Warning: 187 points were rejected as lying outside the specified window
## Warning: data contain duplicated points
marks(bears_all_elevation) <- data.frame(season = marks(bears_all_elevation))
model_poly_elevation <- ppm(bears_all_elevation, ~ polynom(elev, 2) + marks,
covariates = list(elev = elev_clean))
print(model_poly_elevation)
## Nonstationary multitype Poisson process
## Fitted to point pattern dataset 'bears_all_elevation'
##
## Possible marks: 'autumn', 'spring', 'summer' and 'winter'
##
## Log intensity: ~elev + I(elev^2) + marks
##
## Fitted trend coefficients:
## (Intercept) elev I(elev^2) marksspring markssummer
## -1.850394e+01 -3.937116e-03 9.510112e-07 -3.957121e-02 7.891398e-01
## markswinter
## -2.379296e+00
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -1.850394e+01 4.652449e-02 -1.859513e+01 -1.841276e+01 ***
## elev -3.937116e-03 1.034097e-04 -4.139795e-03 -3.734437e-03 ***
## I(elev^2) 9.510112e-07 6.143974e-08 8.305915e-07 1.071431e-06 ***
## marksspring -3.957121e-02 5.053363e-02 -1.386153e-01 5.947288e-02
## markssummer 7.891398e-01 4.266227e-02 7.055233e-01 8.727563e-01 ***
## markswinter -2.379296e+00 1.215116e-01 -2.617454e+00 -2.141137e+00 ***
## Zval
## (Intercept) -397.7247657
## elev -38.0729907
## I(elev^2) 15.4787625
## marksspring -0.7830669
## markssummer 18.4973701
## markswinter -19.5808065
# Compares models using AIC
cat("AIC for poly model: ", AIC(model_poly_elevation), "\n")
## AIC for poly model: 142038.1
cat("AIC for gam model: ", AIC(gam_elevation), "\n")
## AIC for gam model: 4186.001
Model shows a negative relationship between bear locations and elevation, suggesting bears are less likely to be found at higher elevations. Quadratic term suggests a slight non-linear effect, but primary relationship appears negative. Model also indicates that summer sees the highest concentration of bear points, while winter has an expected decrease in bear points. Spring also shows a negative effect, although it is smaller than winter’s effect.
Bears tend to prefer lower elevations (based on the negative relationship with elevation). The summer season has the highest bear presence, while winter and spring have much lower bear occurrences, as to be expected.
# Creates seasonal ppp objects
ppp_list <- list()
for (s in names(seasons)) {
bears_loc_season <- bears_loc_filtered[bears_loc_filtered$month %in% seasons[[s]], ]
bears_sf_season <- st_as_sf(bears_loc_season, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)
bears_sf_season_proj <- st_transform(bears_sf_season, crs = 3005)
coords_season <- st_coordinates(bears_sf_season_proj)
suppressWarnings(
ppp_list[[s]] <- ppp(x = coords_season[,1], y = coords_season[,2], window = window)
)
}
# Extracts covariates
elevation <- DATA$Elevation
cover <- DATA$Forest
hfi <- DATA$HFI
dist_water <- DATA$Dist_Water
# KDE comparison
dens_list <- lapply(ppp_list, function(p) density(p, sigma = 10000))
for (s in names(dens_list)) {
plot(dens_list[[s]],
main = paste("KDE Surface -", s))
}
# Quadrat test per season
for (s in names(ppp_list)) {
qt <- quadrat.test(ppp_list[[s]], nx = 5, ny = 5)
print(qt)
plot(qt,
main = paste("Quadrat test -", s))
}
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of CSR using quadrat counts
##
## data: ppp_list[[s]]
## X2 = 561.39, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 21 tiles (irregular windows)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of CSR using quadrat counts
##
## data: ppp_list[[s]]
## X2 = 2202.7, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 21 tiles (irregular windows)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of CSR using quadrat counts
##
## data: ppp_list[[s]]
## X2 = 3767, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 21 tiles (irregular windows)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of CSR using quadrat counts
##
## data: ppp_list[[s]]
## X2 = 2720.5, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 21 tiles (irregular windows)
Dividing area into a grid of 5x5 and testing for deviations from a random distribution. Bear locations do not follow a random distribution and show spatial clustering across all seasons, meaning their presence is likely influenced by specific environmental or ecological factors that are season-dependent. The quadrant test results suggest that for all seasons (winter, spring, summer, and autumn), bear locations exhibit significant spatial clustering or non-random patterns rather than being randomly distributed (as indicated by the very small p-values).
Bear distribution is influenced by factors other than randomness, such as environmental features, resource availability, or behavior patterns that vary across seasons.
Winter: X2 = 561.39, p-value < 2.2e-16. P-value indicates a strong deviation from CSR for winter, suggesting that bear locations are not randomly distributed in the winter season.
Spring: X2 = 2202.7, p-value < 2.2e-16. Test statistic is very large, and the p-value is highly significant, indicating that the bear locations in spring also show non-random clustering or patterns.
Summer: X2 = 3767, p-value < 2.2e-16. Summer further supports a non-random pattern, with an even higher test statistic and significant p-value.
Autumn: X2 = 2720.5, p-value < 2.2e-16. Autumn shows a significant departure from CSR, confirming that bear locations are non-randomly distributed.
# Inhomogeneous Poisson model per season
for (s in names(ppp_list)) {
model <- ppm(ppp_list[[s]] ~ elev, covariates = list(elev = elevation))
plot(predict(model),
main = paste("Intensity ~ Elevation -", s))
}
## Warning: glm.fit: algorithm did not converge
## Warning: Values of the covariate 'elev' were NA or undefined at 0.02% (1 out of
## 5456) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## ppp_list[[s]], trend = ~elev, data = NULL, interaction = NULL,
# Multivariate ppm
for (s in names(ppp_list)) {
model <- ppm(ppp_list[[s]] ~ elev + cover + hfi + dist_water,
covariates = list(elev = elevation, cover = cover, hfi = hfi, dist_water = dist_water))
plot(predict(model),
main = paste("Predicted Intensity -", s))
}
## Warning: Values of the covariate 'hfi' were NA or undefined at 0.68% (4 out of
## 587) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## ppp_list[[s]], trend = ~elev + cover + hfi + dist_water,
## Warning: glm.fit: algorithm did not converge
## Warning: Values of the covariate 'hfi' were NA or undefined at 0.45% (12 out of
## 2667) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## ppp_list[[s]], trend = ~elev + cover + hfi + dist_water,
## Warning: Values of the covariates 'elev', 'hfi' were NA or undefined at 0.4%
## (22 out of 5456) of the quadrature points. Occurred while executing: ppm.ppp(Q
## = ppp_list[[s]], trend = ~elev + cover + hfi + dist_water,
## Warning: Values of the covariate 'hfi' were NA or undefined at 0.67% (18 out of
## 2698) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## ppp_list[[s]], trend = ~elev + cover + hfi + dist_water,
For each season - assessing spatial distribution of bear occurrences (I.e., clustered, dispersed, or randomly distributed)
# Ripley's K with Envelope for each season
for (s in names(ppp_list)) {
K <- Kest(ppp_list[[s]])
E <- envelope(ppp_list[[s]],
Kest,
correction = "border",
rank = 1,
nsim = 19,
fix.n = TRUE)
plot(K,
main = paste("Ripley's K -", s))
plot(E,
add = TRUE,
lty = 2)
}
## Generating 19 simulations of CSR with fixed number of points ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
## Generating 19 simulations of CSR with fixed number of points ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
## Generating 19 simulations of CSR with fixed number of points ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
## Generating 19 simulations of CSR with fixed number of points ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
Across all four seasons, the Ripley’s K function plots show a
pattern of clustering. This can be concluded as the estimated K
functions (iso, trans, bord) are
generally significantly above the expected K function under complete
spatial randomness and they often lie outside the simulation
envelope.
Winter: Clustering may reflect limited movement due to snow, or congregation around scarce food sources during milder periods. It could also be indirectly related to dens.
Spring: Clustering could indicate post-hibernation dispersal patterns, concentration around newly available food sources, habitat preferences, or potential attraction to human-related food.
Summer: Clustering may reflect the distribution of abundant food resources, proximity to water sources, preferred habitats, and ongoing attraction to human-related food during peak bear activity.
Autumn: Clustering could reflect bears concentrating on late-season food sources for hibernation, movements related to finding den sites, and continued attraction to human-related food as natural food availability changes.